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Abstract 

Supersonic flow of Bose-Einstein condensate past macroscopic obstacles is studied 
theoretically. It is shown that in the case of large obstacles the Cherenkov cone trans- 
forms into a stationary spatial shock wave which consists of a number of spatial dark 
solitons. Analytical theory is developed for the case of obstacles having a form of a 
slender body. This theory explains qualitatively the properties of such shocks observed 
in recent experiments on nonlinear dynamics of condensates of dilute alkali gases. 

In usual compressible fluid dynamics, there are two typical situations when shock waves 
can be generated. The first one is connected with breaking of a nonlinear wave and the second 
with a supersonic flow past a body (see, e.g., [TJ El)- In a viscous fluid, the shock wave 
can be represented as a narrow region within which strong dissipation processes take place 
and the thermodynamic and hydrodynamic parameters of the flow undergo sharp change. 
If viscosity is negligibly small compared with dispersion effects, the shock discontinuity is 
resolved into an expanding region filled with nonlinear oscillations. A remarkable feature 
of such a "dispersive shock" is generation of solitons at one of its boundaries so that the 
whole structure can often be asymptotically represented as a "soliton train". The theory 
of dispersive shocks based on the Whitham nonlinear modulation theory was developed for 
media described by the Korteweg-de Vries (KdV) equation as for the wave breaking case jl] , 
so for the flow past a slender body In the latter case the "soliton train" represents a 
"fan" of oblique spatial solitons spreading downstream from the pointed part of the body. 

After experimental discovery of the Bose-Einstein condensate (BEC) jHl 13 IE], it was 
found that its dynamics is described very well by the Gross-Pitaevskii (GP) equation (see, 
e.g., 0) 

it M = + + 

at 2m 

where ^(r) is the order parameter ("condensate wave function"), U(r) is the potential which 
confines atoms of Bose gas in a trap, and g is an effective coupling constant arising due to 



inter-atomic collisions with the s-wave scattering length a s , 

g = A7ih 2 a s /m, (2) 

m being the atomic mass. The GP equation (0) combines the dispersive and nonlinear effects, 
and the corresponding properties of BEC dynamics have been investigated in a number of 
papers (see for review, e.g., ||). In particular, if g > 0, then existence of dispersive shocks 
can be expected, their theory was developed in [TUl E] and they were observed in a recent 
experiment ^21- Although in the experiment ^2] the shock flow was strongly disturbed by 
vortices arising due to rotation of the condensate, we were informed ^3] about unpublished 
results of experiments on shocks in non-rotating BEC, and these results agree qualitatively 
with the theory ^H] • in experiments ^3 EI] , the shocks were generated by sharp push of 
cylindrical BEC from its axis by a laser beam propagating along the axis. After breaking of a 
cylindrical nonlinear wave the dispersive shock occurs which consists of a train of concentric 
cylindrical solitons. 

Besides the experiments on generation of shocks after wave breaking of BEC's flow, in 
the experiments were performed on BEC's flow past an obstacle after release of BEC from 
the confining potential. The problem of superflow past a body has been studied intensively 
in connection with a problem of critical velocity v c above which superfluidity disappears (see, 
e -g- ; It was found that superfluidity is lost due to generation of vortex rings behind 

an obstacle which gives the estimate of critical velocity 

v c ~(H/dm)Hd/£), (3) 

where £ = H/ ^Jlmn^g is the healing length (i.e., "vortex core size"), d is the size of obstacle 
in transverse direction, and n Q is the density of atoms in the condensate far from the obstacle. 
For large obstacles with d 3> £ this estimate gives the critical velocity much less than the 
sound velocity 

c s = h/V2m£. (4) 

This transition to dissipative behaviour in quantum fluids described by the GP equation 
was confirmed by numerical study [23 H3 EE] where it was found that indeed vortices are 
generated at velocities above the critical one about ~ 0.45c s for d = 10£. 

Since the radius of vortex rings (or distance between vortices in pairs in two dimensions) 
is about the obstacle size d, this mechanism of vortices emission becomes ineffective for d ~ £, 
and for such small bodies ( "impurities" ) the critical velocity arises due to Cherenkov emission 
of sound waves and coincides, hence, with the sound velocity c s ^H]- Obviously, emission 
of waves in a supersonic flow past an obstacle remains effective also for large obstacles 
with d ^> £, but in this case emitted waves are not linear sound waves and, instead of 
a Cherenkov cone, we arrive at the above mentioned dispersive shock consisting of oblique 
spatial solitons. Actually, these shock waves have been observed in the experiment ^J] where 
the wave pattern consists of a series of distinct oblique traces which cannot be attributed to 
a linear Cherenkov radiation implying appearance of a single cone but can be explained by 
generation of supersonic spatial solitons in the flow similar to those predicted by the KdV 
dynamics in [3], [Sj- An easy estimate shows that, after long enough time of expansion, the 
velocity of the flow past a body is much greater than the local sound velocity in BEC near 
this body. Indeed, the flow velocity u in a free expansion has the order of magnitude of the 
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sound velocity at the center of BEC before its release and it is known that the sound speed 
in BEC is proportional to the square root of density. Since for the expansion time t ^> ujJ 1 
(u>± is the radial trap frequency before expansion) the density is proportional to t~ 2 (see, 
e.g., ,2(X ), we find that the ratio of the expansion flow velocity to the local sound speed, i.e. 
the Mach number M, is about 

M ~ u±t ^> 1 (5) 

for t ^> u>~]_ . This corresponds approximately to conditions of the experiment ^3] where 
formation of a "fan" of spatial solitons was observed. Motivated by the results of this 
experiment and above physical argumentation, we shall develop here the theory of hypersonic 
spatial dispersive shocks on the basis of the GP equation (JTJ). 

In accordance with the experiment IT3] , we consider a two-dimensional flow of the 
condensate, so that the condensate wave function ip depends on only two spatial coordinates 
r = (x, y). To simplify the theory, we assume that the characteristic size of the body is much 
less than its distance from the center of the trap, so that incoming flow can be considered 
as uniform with constant density tiq of atoms and constant velocity uo directed parallel to 
x axis. It is convenient to transform Eq. JI} to a "hydrodynamic" form by means of the 
substitution 

V>(r,*) = Vn(r,t)exp U j u(r',t)d/) , (6) 

where n(r, t) is density of atoms in BEC and u(r, t) denotes its velocity field, and to introduce 
dimensionless variables x = x/\r2£, y = y/V2£, t = (c a /2\/2£)t, n = n/n , u = u/c s , where 
numerical constants are introduced for future convenience. As a result of this transformation 
we obtain the system (we omit tildes for convenience of the notation) 

\n t + V(rm) = 0, 

i t t~i t~i T(Vn) 2 An 

§11* + (uV)u + Vn + V 



(7) 

V ; 



in 2 An 

(where V = (d x ,d y )) which should be solved with the boundary conditions 

n = 1, u = (M, 0) at x — > — oo (8) 

for incoming flow and 

u-N| 5 = (9) 

at the body surface S, where N denotes a unit vector of outer normal to the surface S. 
Under our assumption that the size of the body is much less than the distance from the 
center of the cylindrically symmetrical flow, the characteristics of the flow near the body 
change with the time very slowly, so that the arising wave pattern can be considered as quasi- 
stationary. Hence, we can confine ourselves to the steady solutions of our problem (JTj)-® 
and replace Eqs. (J7J) by their stationary versions for two-dimensional stationary velocity field 
u = (u(x,y),v(x,y)): 



uu x + vu v + n x + 



x i y I o 2 



nu) x + {nv\, = 0, 

y ^xx H~ T^yy 



n x + 7i n xx + n. 



n 2 An J (10) 



UVx + VVy + ny+ t__JL =0 . 
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If the body is symmetric with respect to x axis and the form of its boundary is given by 
the function y = ±F(x), F(0) = 0, F(L) = 0, L being the longitudinal size of the body in 
dimensionless units, then we can confine ourselves to consideration of only the upper half- 
plane y > 0, so that N oc (F'(x), —1), and the boundary conditions (JSj), are transformed 
to 

n — 1, u — M, v — at x — > — oo, (11) 



v = uF'(x) at y = F(x). (12) 

The system (jlU|) - (|12|) is still too complicated for analytical treatment. However, the flow is 
supposed to be supersonic (see (J5J)) which allows us to asymptotically transform Eqs. (fTUj)- 
(|12|) to a much simpler problem of ID "unsteady" flow along y axis with the scaled x 
coordinate playing the role of "time" |21j. To this end, we substitute into Eqs. (fT0j) new 
variables 

u = M + Ul + 0(1/M), T = x/2M, Y = y, (13) 
assuming M _1 <C 1. Then to leading order we obtain 



\v T + vv Y + n Y + ( - ^ ) =0, 



|n T + (nt> ) Y = 0, 

2 

In 2 An 



(14) 



|ttlT + VUiy = 0. (15) 

Equations ()14j) represent the hydrodynamic form of ID defocusing nonlinear Schrodinger 
(NLS) equation 

i^T + ^YY -2\^\ 2 ^ = (16) 



for a complex field variable 

-Y 



^ = v / ^exp^y v{Y\t)dY^j 



(17) 



and we can apply well-known methods to their study. If n and v are found from the system 
()14j) . then the correction u\ to the a; component of velocity can be calculated with the use 
of Eq. (|15|). It is remarkable that in the case of a slender body, for which Ma < 1 where 
a = max\F'(x)\, the boundary condition (fT2"j) reduces (to leading order in M _1 ) to 

v = Vp = \df/dT at Y = f(T), (18) 

where f(T) = F(2MT), and it does not contain the u- variable in this approximation. 

Thus, we have reduced the problem of flow past a slender body to the classical "piston" 
problem for ID flow along a tube with a piston moving according to the law (|18j) . In ordinary 
gas dynamics, such a piston motion leads to the generation of two shock waves which form 
due to the breaking of the flow profile during two different phases of the piston motion: 
forward and reverse. In the original 2D problem this corresponds to the supersonic flow past 
the front and the rear parts of the body respectively and is accompanied by the occurrence 
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of two spatial shocks (oblique compression jumps) spreading from the body edges (see [I] 
for instance). 

In contrast to the classical gas dynamics, the piston problem is now posed for dispersive 
equations (jl4j) . Before we proceed with the analysis of this problem we outline the qualitative 
structure of the dispersive flow past finite-length body using the results in [3], j^j. This will 
help us then to reformulate the piston problem for the NLS equation in terms of much better 
explored initial-value problem. 

In dispersive hydrodynamics, both shocks (breaking singularities) spreading from the 
body edges resolve into nonlinear oscillatory zones, the dispersive shocks. At finite distances 
from the body surface these two dispersive shocks have similar structure (see each 
represents a modulated nonlinear wave having a form close to a chain of solitons at one edge 
of the oscillatory zone and degenerating into a linear wave at the opposite edge. However, 
at large distances from the body the two dispersive shocks demonstrate drastically different 
behaviour: in the physical systems with negative dispersion studied in |H] the dispersive 
shock spreading from the front edge of the body transforms into a soliton train while the 
dispersive shock forming at the rear end of the body completely degenerates into a small- 
amplitude linear radiation. In the case of the NLS equation describing nonlinear waves with 
positive dispersion the qualitative picture will be reversed, i.e. the "nontrivial" dispersive 
shock transforming into a (dark) soliton train will form due to the flow past the rear part of 
the body (reverse motion of the piston in the corresponding nonstationary problem). In the 




Figure 1: (Y, T)-plane of the piston problem (Y > 0) and the equivalent initial data A+(Y, 0)) 
(Y < 0) for the Riemann invariant. The shaded area marks a "soliton" part of the initial 
profile. Dashed line: the "piston" trajectory Y = f(T). The lines Y* 2 (T), Yi^{T) are the 
edges of the front and rear dispersive shocks respectively. 

present work we will be concerned with this wave only. 

We now formulate the problem of the rear dispersive shock description in a mathemati- 
cally consistent way. The relevant part of the (Y, T)-plane in the auxiliary piston problem 
is subdivided into three distinct regions (see Fig. 1). Generally, for Y > f(T), the "gas" is 
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put into motion by the "piston" moving according to Eq. (|18|) and in the region I near the 
"piston" the gas motion can be described by the dispersionless limit 



\ut + (nv)y — 0, \vt + v vt + ny = (19) 

of Eqs. (jHj) . But formal solution of the dispersionless equations cannot be extended to the 
whole (Y, T)-plane because the F-derivatives blow up along some line in this plane. Hence, 
here we have to take into account the dispersion effects which lead to formation of the region 
// filled with nonlinear oscillations — a dispersive shock. Close to its boundary Y = Yi(T) 
the oscillations tend to a train of dark solitons of the NLS equation and at the opposite 
boundary Y = Y2(T) the amplitude of oscillations tends to zero which corresponds to a linear 
"sound" wave propagating into the undisturbed region III with Y > 3*2 (T). For T>L/Af 
the whole structure can be asymptotically represented as a soliton train Thus, in the 
case of a macroscopic obstacle with characteristic size much greater than the healing length 
£ the Cherenkov cone "unfolds" into a "fan" of solitons. 

The outlined nonlinear dispersive flow is described most conveniently in terms of the 
Riemann invariants (see, e.g., |2"3*1 I22|). In the regions / and /// of the smooth flow these 
are the Riemann invariants of the dispersionless system (|19|) 



which satisfy the equations 



8T = v " r ' ~' 8Y 
where the characteristic velocities are 



\± = \v±^i, (20) 



yA± +y ± (A +) A_)^ = > (21) 



V + = 3A+ + A_ , V^= 3A_ + A + . (22) 

In the dispersive shock region //, there are four Riemann invariants Aj, i = 1,2,3,4, which 
parameterize the modulated periodic solution of the full system (|Hj) . and obey the corre- 
sponding Whitham modulation equations [21] 

<9A <9A 

^+^(Ai,A 2 ,A3,A 4 )^ = 0, i = 1,2,3,4. (23) 

For the case of the defocusing NLS equation the characteristic velocities Vi{\\, . . . , A 4 ) rep- 
resent certain combinations of the complete elliptic integrals of the first and the second kind 
(see also |2*3*1 122j). We do not need the concrete expressions for them here. 

The dispersionless Riemann invariants (J2U|) are constant along the corresponding families 
of characteristics of the system ()19|) and they must satisfy the matching conditions with two of 
the Riemann invariants Aj of the Whitham equations ()23|) along the lines Y = Yi^{T), which 
represent (free) boundaries of the dispersive shock (see 12Z| for a detailed description of 
the matching problem for the solutions of systems (|21|) and ()23|)). 

In the undisturbed region III, where v — 0, n — 1, both Riemann invariants (|2U|) are 
constant: A± = ±1. According to the well-known argumentation about the flow adjacent to 
a simple wave (see pQ, Section 104), one of the Riemann invariants Aj, i = 1,2,3,4, which 
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matches with, say, A_ (for sake of definiteness we assume that it is A 4 ) must also be constant 
in the whole region II: A4 = A_ = —1 in II. On the other hand, the gas motion in the 
region I produced by the "piston" is described by a simple wave solution of Eqs. (fTU|) (see [T] , 
Section 101) again with one Riemann invariant constant everywhere in /. It must match with 
A4 along the characteristic line Y = Yx(T) so that A_ = A4 = —1 in the whole (Y, T)-plane 
including the trajectory of the "piston". Hence, we have at the "piston" v p /2 — ^Jn^, = — 1 
which yields the gas density, n p = (t> p + 2) 2 /4, and, consequently, the values of both Riemann 
invariants: 

A^ = -1, A^ = \df/dT + 1 at Y = f(T). (24) 

To use the method of Ref. (221, we have to transform these boundary conditions at the 
"piston" to the equivalent initial conditions at T = 0, Y < (see Fig. 1). This problem for 
the system (fT$|) can be easily solved using characteristics. Indeed, we have A_ = —1, hence 
A + obeys the simple wave equation following from (J2T|) (see, e.g., (221) 

|± + (3A + - 1)|± = (25) 

whose general solution is 

Y = (3A+-1)T + F(A+). (26) 

The function Y(X + ) must be chosen so that the condition (|2^j) is satisfied. This gives at 
once 

F(A+) = f(r) - (3A + - l)r, (27) 

where r is determined implicitly by the equation A + = \f'{r) + 1. Thus, we arrive at a 
parametric form of the equivalent initial distribution of the Riemann invariant A + : 

A+ = §/'(r) + 1, Y = f(r) - (§/'(r) + 2)r. (28) 

In principle, knowing the initial data for the Riemann invariants X±, one can construct the 
full solution of the Whitham system f)23j) using the extended Gurevich-Pitaevskii problem 
formulation (see [2B],[22]) and thus, to asymptotically describe the dispersive shock for all 
Y, T. However, if one is interested in the asymptotic structure of the flow in the region 
far enough from the body where the rear dispersive shock develops into a "fan" of spatial 
solitons well separated from each other, one can take advantage of a more simple asymptotic 
method of Ref . [221 - 

According to the asymptotic method developed in Ref. |22| each soliton in the soliton 
train evolving from the initial pulse is parameterized by the eigenvalue A& of the generalized 
Bohr-Sommerfeld quantization rule 

j - A + )(A fe - A_) dY = 2n(k + i), k = Q,l,...,K, (29) 

where in our case \+{Y) is given by Eq. (|28J). A_ = —1, and the integration is taken over 
the cycle around two turning points where the integrand functions vanishes. Returning to 
the spatial coordinates (fT3~|). we find the profile of the Afc-soliton in the train as (see [22]) 

1 — A 2 

n k (x, y) = l „ . h , (30) 

cosh 2 [yr^Af( 2/ -(A fc /M)x)]' 
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that is the "fan" of spatial dark solitons in the shock is made of soliton "feathers" lying 
asymptotically along the lines 

y=(\ k /M)x, k = 0,l,...,K, (31) 

in the upper half-plane and symmetric "fan" of solitons is generated in the lower half-plane. 
Let us illustrate this theory by an example of a body with a parabolic profile 

y = F(x) = ax(L -x)/(2M) 2 , < x < L, (32) 

or 

Y = aT(T — T), < T < T , (33) 
where T = L/2M, so that the initial condition (|28j) takes the form 

A + = \a{T Q - 2r) + 1, Y = r(2ar — aT /2 — 2), (34) 

for < t < To, that is in the interval — Tq(2 — 3aTo/2) < Y < 0, and A + = 1 outside this 
interval; see Fig. 2. Its part with A + > 1 evolves into non-solitonic wave which does not give 



J 1 y 

Y 2r 



Figure 2: The plot of A + as a function of Y for the parabolic body profile ([32)1 . Dashed lines 
labeled by A , Ai, A 2 etc denote the levels corresponding to these eigenvalues, Y 2 i and Y 2r 
illustrate positions of the left and right turning points corresponding to the eigenvalue A 2 in 
this instance. 

any contribution into the shock. However, its "well" part A + < 1 leads to the bound states 
in the spectral problem (|2~9j) and, hence, to the train of spatial solitons (jSUj) generated in the 
shock. Integral in ()29|) is calculated in a closed form which gives the equation 

+ 4A fc - f «T )(i«T - 1 + A fc ) 3/2 = k + l k = 0,l,...,K, (35) 

and its roots A& must lie in the interval 1 — aT /2 < A& < 1. The greatest root A^ has 
a value close to unity so that the number of solitons in the shock is given approximately 
by (|33|) with Afc = 1, k — K. To transform this estimate to dimensional variables, we take 
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I = y/2£L = 2V2£MT longitudinal size of the obstacle, d = 2y/2£f(T /2) = \/2£<*T 2 
as its transverse size and obtain 



K = — \l 

3n 



27 Md 

To~T 



id 



(36) 



Although this formula is derived for a slender body, we can use it as a rough estimate of a 
number of solitons in the shock generated in the supersonic flow past an obstacle: 



K = const ■ {Id/ 



2x1/2 



(37) 



where K must be much greater than unity. The most shallow dark soliton lies near the 
Cherenkov cone with angle 9 = 1/M = c s /uq with respect to the direction of the flow. The 
resulting pattern is shown in Fig. 3. 



100 




-100 



Figure 3: The pattern of dispersive shocks generated by a BEC supersonic flow with M = 5 
past a slender body (black cigar-shaped figure). The oblique lines represent the traces of 
dark solitons in the (x, |/)-plane. All dimensions are measured in units of \/2£. 

In conclusion, we have developed the theory of spatial dispersive shock waves generated 
by a flow of Bose-Einstein condensate past a slender body. The theory predicts formation of 
a series of oblique spatial solitons in the flow and explains qualitatively the shock patterns 
observed in the experiment il3J. 

We are grateful to E.A. Cornell and P. Engels for information about the results of the 
experiments of JILA group prior to publication. AMK thanks RFBR (grant 05-02-17351) 
for financial support. 
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